4. Data Exploration

1 Why data exploration and descriptive statistics?

1. Reporting on results: not useful to only report a table with the results for each subject
2. Let the data speak by summarizing and visualizing it
3. Get insight in the data
4. Discover errors, anomalies or even fraud
5. Check assumptions that are required for the downstream statistical inference, e.g. are the data normally distributed

2 Univariate exploration of quantitative variables

2.1 Histogram

NHANES %>% filter(Gender=="female") %>% 
ggplot(aes(x=DirectChol)) + 
geom_histogram(aes(y=..density.., fill=..count..),bins=30) +
geom_density(aes(y=..density..))

1. Select females and pipe results to ggplot

   NHANES %>% filter(Gender=="female")

2. Select data to plot

ggplot(aes(x=DirectChol)) +

3. Equal bins for interpretation, number of bins can be selected with the bins argument to the geom_hist

4. Relative frequenties to enable visual comparison between histograms

geom_histogram(aes(y=..density.., fill=..count..)) +

5. If we have enough observations we can use a kernel density estimator of f(x)

geom_density(aes(y=..density..))

2.2 Boxplot

With ggplot we always have to define an x variable if we make a boxplot. If we use a string then all data is put in one category. And one boxplot is constructed.

NHANES %>% filter(Gender=="female") %>%
  ggplot(aes(x="",y=DirectChol)) +
  geom_boxplot() 

So we can add a boxplot to a ggplot figure by using the geom_boxplot function.

If the dataset is small to moderate in size we can also add the raw data to the plot with the geom_point() function and the position="jitter" argument. Note, that we then also set the outlier.shape argument in the geom_boxplot function on NA so that the outliers are not plotted twice.

Here, we will plot again the relative abundances of Staphylococcus from the armpit transplant experiment

ap<-read_csv("https://raw.githubusercontent.com/GTPB/PSLS20/master/data/armpit.csv")
ap

# A tibble: 20 x 2
   trt          rel
   <chr>      <dbl>
 1 placebo     55.0
 2 placebo     31.8
 3 placebo     41.1
 4 placebo     59.5
 5 placebo     63.6
 6 placebo     41.5
 7 placebo     30.4
 8 placebo     43.0
 9 placebo     41.7
10 placebo     33.9
11 transplant  57.2
12 transplant  72.5
13 transplant  61.9
14 transplant  56.7
15 transplant  76  
16 transplant  71.7
17 transplant  57.8
18 transplant  65.1
19 transplant  67.5
20 transplant  77.6

ap %>%  ggplot(aes(x=trt,y=rel)) + geom_boxplot(outlier.shape=NA) + geom_point(position="jitter")

When we specify a factor variable for x, we get a boxplot for each treatment group.

2.3 Descriptive statistics

2.3.1 Central location: Mean or Median?

2.3.1.1 Mean

• In a period of 30 years males hope to have on average 64.3 partners and females 2.8

2.3.1.2 Median

• The median of the number of partners males and females want to have is both 1

2.3.1.3 What happens?

desired number of partners

• Mean is very sensitive towards outliers!

2.3.2 Geometric mean

\[\sqrt[n]{\prod\limits_{i=1}^n x_i} = \exp\left\{\frac{1}{n} \sum_{i=1}^n \log(x_i)\right\}\]

• Geometric mean is closer to the mediaan then the mean

• log-transformation removes skewness

• Often a more useful measure for the central location than median:

1. Uses all observations: is more precise
2. It is the ordinary mean on log-transformed data \(\rightarrow\) classical statistical methods can be directly applied, e.g. hypothesis tests and confidence intervals (see chapter 5)
3. Useful for many biological characteristics e.g. concentrations that cannot be negative.
4. Differences on a log scale have the interpretation of a log fold change:

\[\log (B) - \log(A)= \log(\frac{B}{A})=\log(FC_\text{B vs A})\]

logSummary <- 
NHANES %>% filter(Gender=="female") %>% summarize(logMean=mean(DirectChol %>% log,na.rm=TRUE),sd=sd(DirectChol %>% log,na.rm=TRUE),mean=mean(DirectChol,na.rm=TRUE),median=median(DirectChol,na.rm=TRUE)) %>% mutate(geoMean=exp(logMean))

NHANES %>% filter(Gender=="female") %>%
ggplot(aes(x=DirectChol %>% log)) + 
geom_histogram(aes(y=..density.., fill=..count..),bins=30) +
geom_density(aes(y=..density..)) +
  stat_function(fun=dnorm,color="red",args=list(mean=logSummary$mean, sd=logSummary$sd))

logSummary

    logMean        sd     mean median  geoMean
1 0.3564553 0.2731811 1.482467   1.42 1.428258

• Indeed the mean is pulled to larger values by the skewed data
• The geometric mean is closer to the median
• The cholestorol data are much more symmetric upon log transformation and the approximation by a normal distribution is good.

2.4 Descriptives for the variability

The variability around the central value is crucial:

1. Biologists are often interested in how animals or plants are spread in the study region.
2. Compare groups: the group effect is more clear when the response has less variability. Quantifying variability is crucial to distinguish between systematic and random patterns.

• The response varies between and within individuals and is the reason why we need statistics.

• Crucial to describe both the central location and the variability.

• Which part of the variability can we explain (e.g. with characteristics treatment, age, etc,) and which part is unexplained? ### Sample variance and sample standard deviation

• Sample variance: \[ s_X^2= \sum\limits_{i=1}^n \frac{(X-\bar X)^2}{n-1} \]

• Interpretation is often difficult because it is in another unit than the measurements

• Standard deviation: \[ s_X= \sqrt{s_x^2} \]

• Very useful for Normal distributed observations:

  • 68% of the observations falls in the interval \(\bar{x} - s_x\) en \(\bar{x} + s_x\)
  • 95% of the observations falls in the interval \(\bar{x} - 2 s_x\) en \(\bar{x} + 2 s_x\).

• These intervals are referred to as 68% en 95% reference-intervals.

• If the data are not normally distributed reference intervals are not valid

2.4.1 Interquartile range

For skewed data the standard deviation is not useful - It is very sensitive to outliers - Inter Quartile Range: Distance between first and third quartile - Width of the boxplot!

NHANES %>% filter(Gender=="female") %>% summarize(IQR=IQR(DirectChol,na.rm=TRUE))

   IQR
1 0.52

NHANES %>% filter(Gender=="female") %>%
  ggplot(aes(x="",y=DirectChol)) +
  geom_boxplot() 

3 Normale approximation

• Biological and chemical data are often normally distributed upon transformation.

• If this is the case we can get a lot of insight in the data using just to descriptive statistics: gemiddelde \(\mu\) en standaard deviatie \(\sigma\).

3.1 Evaluation with QQ-plots

If your analysis build upon the assumption that the data are Normally distributed it has to be verified.

We use QQ-plots or quantile-quantile plots).

• Observed quantiles from the observations in the sample are plotted against quantiles from the Normale distribution.

• If the data are normally distributed both quantiles have to be in line.

• Dots in the plot are expected on a straight line

• Systematic deviations of the straight line indicate that the data are not normally distributed.

• Note, that we will always observe some random deviations from the straight line in the plot because of random biological variability, which is not indicative for deviations from Normality.

• So it is important to train yourself to learn to distinguish systematic from random deviations.

3.1.1 Normal data

• We will first simulate data from the normal distribution to show how the plots look like for data that is meeting the assumptions.

• We will simulate data from 9 samples with a mean of 18 and standard deviation of 9

n <- 20
mu <- 18
sigma <- 9
nSamp <-9

normSim <- matrix(rnorm(n*nSamp,mean=mu,sd=sigma),nrow=n) %>% as.data.frame

normSim %>% gather(samp,data) %>% 
ggplot(aes(x=data)) +
geom_histogram(aes(y=..density.., fill=..count..),bins=30) +
geom_density(aes(y=..density..)) +
  facet_wrap(~samp)

normSim %>% gather(samp,data) %>% 
  ggplot(aes(sample=data)) +
  geom_qq() +
  geom_qq_line() +
facet_wrap(~samp)

So even for normal data we observe some deviations!

3.1.2 Real data

NHANES %>% filter(Gender=="female"&!is.na(BMI)) %>%
  ggplot(aes(x=BMI))+
   geom_histogram(aes(y=..density.., fill=..count..)) +
   xlab("BMI") +
   ggtitle("All females in study") +
  geom_density(aes(y=..density..))

NHANES %>% filter(Gender=="female"&!is.na(BMI)) %>%
  ggplot(aes(sample=BMI)) +
  geom_qq() +
  geom_qq_line()

The QQ-plot shows that the quantiles of the data - are larger (above the line) than these from the Normal in the left tail: compression of the lower tail - are larger (above the line) than these from the Normal in the right tail: long tail to the right - We can clearly see that the data are skewed